GENERALIDADES
La regresión logística es un tipo de análisis de regresión utilizado para predecir el resultado de una variable dependiente o respuesta categórica (por lo general de dos categorías o binaria) en función de una o más variables independientes o predictoras.
La regresión logística es uno de los procedimientos estadísticos más usados en ciencias médicas, epidemiología y ciencias sociales, entre otras.
HISTORIA
La técnica de la regresión logística se originó en el año 1961 con el trabajo de Cornfield, Gordon y Smit
Su uso se incrementa desde principios de los 80 como consecuencia de los adelantos ocurridos en el campo de la computació
¿Por qué regresión logística y no lineal?
Si una variable cualitativa con dos niveles se codifica como 1 y 0, matemáticamente es posible ajustar un modelo de regresión lineal por mínimos cuadrados \(\beta_0+\beta_1X\) (a+bx).
El problema de esta aproximación es que, al tratarse de una recta, para valores extremos del predictor, se obtienen valores de Y menores que 0 o mayores que 1, lo que entra en contradicción con el hecho de que las probabilidades siempre están dentro del rango [0,1]
Para evitar estos problemas, la regresión logística transforma el valor devuelto por la regresión lineal \(\beta_0+\beta_1X\) empleando una función cuyo resultado está siempre comprendido entre 0 y 1.
Existen varias funciones que cumplen esta descripción, una de las más utilizadas es la función logística (también conocida como función sigmoide):
\[\huge P_(Y=k|X=x)=\frac{e^{\beta_0+\beta_1X}}{1+e^{\beta_0+\beta_1X}}\]
Concepto de ODDS o razón de probabilidad
En la regresión lineal simple, se modela el valor de la variable dependiente Y en función del valor de la variable independiente X.
Sin embargo, en la regresión logística, se modela la probabilidad de que la variable respuesta Y pertenezca al nivel de referencia en función del valor que adquieran los predictores, mediante el uso de LOG of ODDs.
Supóngase un evento cuya probabilidad de ser verdadero es de 0.8 (p=0,8). Dado que es una variable binaria, la probabilidad ser falso es de 1 - 0.8 = 0.2 (q=0.2).
Los odds o razón de probabilidad se definen como el ratio entre la probabilidad de evento verdadero y la probabilidad de evento falso p/q. En este caso los odds son 0.8 / 0.2 = 4, lo que equivale a decir que se esperan 4 eventos verdaderos por cada evento falso.
Coeficientes de los predictores
La estimación de los parámetros β0 y β1 se obtiene por el método de máxima verosimilitud (maximum likelihood ML), es decir, el valor de los parámetros β0 y β1 con los que se maximiza la probabilidad de obtener los datos observados.
El método de maximum likelihood está ampliamente extendido en la estadística pero su implementación y explicación no siempre es fácil de realizar y en nuestra caso particular, excede el nivel de este curso.
Afortunadamente, programas como SPSS y R, entre otros, calculan los coeficientes de forma automática.
Interpretación del modelo
Los principales elementos que hay que interpretar en un modelo de regresión logística son los coeficientes de los predictores:
β0 = es la ordenada en el origen o intercept. Se corresponde con el valor esperado del logaritmo de odds cuando todos los predictores son cero. Puede transformarse a probabilidad con la ecuación \(𝑒^{𝛽_0}/(1+𝑒^{𝛽_0} )\). Tras la transformación, su valor se corresponde con la probabilidad esperada de pertenecer a la clase 1 cuando todos los predictores son cero.
βi = los coeficientes de regresión parcial de cada predictor indican el cambio promedio del logaritmo de odds al incrementar en una unidad la variable predictora Xi, manteniéndose constantes el resto de variables. Esto equivale a decir que, por cada unidad que se incrementa Xi, se multiplican los odds por \(𝑒^{𝛽_i}\)
Condiciones para su Implementación
Independencia: las observaciones tienen que ser independientes unas de otras.
Relación lineal entre el logaritmo natural de odds y la variable continua: patrones en forma de “U” son una clara violación de esta condición.
La regresión logística no precisa de una distribución normal de la variable continua independiente.
Número de observaciones: no existe una norma establecida al respecto, pero se recomienda entre 50 a 100 observaciones. Por su parte, Freeman en 198
Predicciones
Una vez estimados los coeficientes del modelo logístico, es posible conocer la probabilidad de que la variable dependiente pertenezca al nivel de referencia, dado un determinado valor del predictor. Para ello se emplea la ecuación del modelo:
\[\Large P_{(Y=1|X)}=\frac{e^{\beta_0+\beta_1X}}{1+e^{\beta_0+\beta_1X}}\] Convertir probabilidad en clasificación
Una de las principales aplicaciones de un modelo de regresión logística es clasificar la variable cualitativa en función de valor que tome el predictor.
Para conseguir esta clasificación, es necesario establecer un threshold (umbral) de probabilidad a partir de la cual se considera que la variable pertenece a uno de los niveles. Por ejemplo, se puede asignar una observación al grupo 1 si \(\hat p_{(Y=1|X)}>0.5\) y al grupo 0 si es lo contrario.
EJEMPLO
El fichero craneos.xlsx contiene los datos recolectados por el Coronel L.A. Waddel en el sureste y este del Tíbet y reportando por G.M. Morant en 192
Se cree posible que los dos conjuntos pertenezcan a dos etnias diferentes. Cada uno de los 32 craneos tienen 5 medidas Longitud (logitud mayor del cráneo), Anchura (anchura mayor horizontal del cráneo), Altura (altura del cráneo), Altura_cara (altura de la cara superior), Anchura_cara (ancho del cara, entre los puntos más externos de los pómulos)
El link de los datos originales los puede encontrar aqui: https://math.univ-angers.fr/~labatte/enseignement%20UFR/master%20MIM/skull.pdf
Cargamos los datos
library(readxl) #Cargamos la librería
datos <- read_excel("data/craneos.xlsx")
str(datos) #Revisamos el tipo de datos
tibble [32 x 6] (S3: tbl_df/tbl/data.frame) $ Longitud : num [1:32] 19050 17250 16700 16950 17500 ... $ Anchura : num [1:32] 15250 13200 13000 15050 13850 ... $ Altura : num [1:32] 14500 12550 12550 13350 12600 ... $ Altura_Cara : num [1:32] 7350 6300 6950 6450 7750 7150 7050 7350 7000 6200 ... $ Anchura_Cara: num [1:32] 13650 12100 11950 12800 13550 ... $ Tipo_raza : num [1:32] 1 1 1 1 1 1 1 1 1 1 ...
En seguida, tenemos que transformar la variable dependiente a 0 y 1 (asi lo requieren las regresiones logisticas binomiales)
datos$Tipo_raza <- ifelse(datos$Tipo_raza==1,0,1) table(datos$Tipo_raza) #verificamos el cambio
0 1 17 15
Para la realización de un modelo logístico, usaremos la función glm(
logist <- glm( Tipo_raza ~ Altura_Cara+Anchura_Cara, data = datos, family = binomial) summary(logist)
Obtener los Odds ratio
Para conocer el impacto que tiene una variable independiente sobre una de la categorías de la variable dependiente en comparación con el otro, usaremos los odds ratio.
Para calcular los odds ratio en r se usa el siguiente comando
exp(coef(x))
En donde x es el objeto que contiene el modelo de la regresión logística, y exp(), calcula el valor exponencial de un número. Esto es equivlente a \(e^x\), donde e es la constante de Euler = 2.71828…
Por ejemplo:
exp(coef(logist))
(Intercept) Altura_Cara Anchura_Cara 3.030043e-16 1.003117e+00 1.000951e+00
Estos Odd, representan una medida relativa de cambio en el riesgo de un perfil respecto del otro.Para interpretarlos, consideraremos los siguiente:
En este caso particular, podemos interpretar que la posibilidad que un cráneo pertenezca a la etnia Tipo 2 (codigo 1 en nuestra base), aumenta 1.0031 veces cuando se aumenta la altura de la cara en una unidad (1 mm).
Tambien podemos interpretarlo como probabilidad: (1.0031 - 1)*100 = 0.31%
paste(round((exp(coef(logist))[2]-1)*100,2), "%")
[1] "0.31 %"
Es decir, la posibilidad de pertenecer a la etnia Tipo 2 aumenta en 0.31% cuando la altura de la cara aumenta en una unidad (1 mm).
Ahora bien, si aumentamos en 10 mm (1 cms) la altura de la cara, la posibilidad de pertenecer a la etnia Tipo 2 se incrementa en: (exp(0.003112088 * 10) -1 )*100 = 3.16%
(exp(logist$coefficients[2]*10) -1)*100
Obtener los Odds ratio
Para la Anchura de la cara, hacemos el mismo procedimiento, la posibilidad de que un cráneo pertenezca a la etnia Tipo 2, aumenta 1.0009 veces cuando se aumenta la anchura de la cara en una unidad (1 mm).
Tambien podemos interpretarlo como probabilidad: (1.00095 - 1)*100 = 0.095%
print(paste((exp(coef(logist))[3]-1)*100, "%"))
[1] "0.0951171037421927 %"
En este caso, vemos que la contrubución del ancho de cara es muy poca y por algo, no es una variable significativa para el modelo.
En conclusión, Altura_Cara aporta significativamente en la estimación del Tipo o afinidad racial, pero Anchura_Cara no.
Cuando los odds ratio es:
Como ya conocemos los coeficientes, podemos construir la ecuación regresiva y calcular posteriormente la estimación de la clasificación.
Sabemos que la ecuación de regresión es: \(\hat Y= \beta_0+\beta_1X_1+\beta_2X_2\), por tanto sustituyendo en la formula, quedaría asi: \[ \hat Y = -35.37278+(0.00311*Altura Cara)+(0.00095*Anchura Cara)\]
La clasificación estimada la podemos realizar con la siguiente función:
predictor <- predict(logist, newdata = datos, type = "response") #Esta función solo nos produce la probabildad de pertenencia head(as.data.frame(predictor))
predictor 1 0.529559315 2 0.009728385 3 0.060501271 4 0.029580447 5 0.780423933 6 0.263678223
Ahora tenemos que clasificar según el umbral que nosotros, decidamos, en este caso 0.5
clasificador <- ifelse(predictor > 0.5, 1, 0) head(as.data.frame(clasificador))
clasificador 1 1 2 0 3 0 4 0 5 1 6 0
Ahora podemos tabular la clasificación estimada con la clasificación observada y ver las coincidencias
matriz_confusion <- table(observaciones = datos$Tipo_raza, Prediccion = clasificador) matriz_confusion
Prediccion
observaciones 0 1
0 13 4
1 4 11
Luego calculamos la Exactitud (Accuracy) como: (13+11/32)*100 = 75%. En R podemos usar este comando:
mean(clasificador==datos$Tipo_raza)*100
[1] 75
En resumen, el modelo es capaz de clasificar correctamente el 75% de las observaciones
Finalmente realizamos la evaluación del modelo con el Test Pseudo R2 de McFadden. Es una medida de bondad de ajuste de un modelo. Se define como la proporción de variación de la variable respuesta o dependiente que es explicado por el modelo de regresión ajustado.
library(DescTools) PseudoR2(logist)
McFadden 0.3703261
Hay varias formas de graficas, pero la más usual es usar ggplot. En nuestro caso, como solo Altura_Cara, fué significativa, procedemos a trabajar son esa variable, ya que con dos o mas variables, es más complicado realizar el gráfico.
library(ggplot2)
ggplot(datos, aes(x=Altura_Cara, y=Tipo_raza)) + geom_point(aes(color = as.factor(Tipo_raza)),
shape = 16, size=5)+
geom_smooth(method = "glm", method.args = list(family = "binomial"), color = "gray20", se = FALSE) +
theme_bw()+ theme(legend.position = "none")
El Análisis Discriminante Lineal o Linear Discrimiant Analysis (LDA) es un método de clasificación supervisado de variables cualitativas en el que dos o más grupos son conocidos a priori y nuevas observaciones se clasifican en uno de ellos en función de sus características. Haciendo uso del teorema de Bayes, LDA estima la probabilidad de que una observación, dado un determinado valor de los predictores, pertenezca a cada una de las clases de la variable cualitativa, P(Y=k|X=x). Finalmente se asigna la observación a la clase k para la que la probabilidad predicha es mayor.
Es una alternativa a la regresión logística cuando la variable cualitativa tiene más de dos niveles.
Cuando se trata de un problema de clasificación con solo dos niveles, ambos métodos suelen llegar a resultados similares.
Ha sido muy útil en biología para poder conocer un set de variables morfológicas que sean capaces de discriminar por ejemplo especies, dimorfismo sexual, etc.
En antropología forense han sido muy utilizadas para poder conocer el sexo o la ancestría de un individuo a partir de su morfología esqueletal.
Esencialmente, los análisis discriminantes buscan la mejor combinación de variables independientes para maximizar la distancia entre grupos.
El proceso de un análisis discriminante puede resumirse en 6 pasos:
El análisis discriminante lineal se utiliza con frecuencia como técnica de reducción de dimensionalidad para el reconocimiento o clasificación de patrones.
Los objetivos del Análisis Discriminante pueden sintetizarse en dos:
La base del análisis discriminante consiste en llegar a establecer una combinación lineal de las variables independientes que nos permita clasificar a los sujetos (poblaciones) en los diferentes grupos establecidos a priori
Esta combinación lineal es lo que se conoce como función discriminante
Hay tantas funciones discriminantes como grupos menos 1
El origen del análisis discriminante lineal se remonta al matemático R. Fisher (1936
Fisher planteó una aproximación en la que el espacio p-dimensional (donde p es el número de predictores originales) se reduce a un subespacio de menos dimensiones formado por las combinaciones lineales de las variables originales que mejor explican la separación de las clases. Una vez encontradas dichas combinaciones, se realiza la clasificación en este subespacio. El subespacio óptimo, según Fisher, es aquel que maximiza la distancia entre grupos en términos de varianza.
Los análisis discriminantes buscan principalmente:
Encontrar relaciones lineales entre las variables continuas que mejor discriminen en los grupos dados a las observaciones.
Construir una regla de decisión que asigne una observción nueva con un cierto grado de riesgo, cuya clasificación previa se desconoce, a uno de los grupos prefijados (como en el caso de la antropología forense)
Limitaciones
Están diseñados para variables independientes que deben distribuir normalmente, lo que lo convierte en una estrategia más limitada que las regresiones logísticas
No solo exigen normalidad de cada variable independiente sino también normalidad multivariante de todas las variables independientes involucradas. Esto convierte a los análisis de funciones discriminantes en una técnica muy limitada ya que es muy difícil encontrar normalidad multivariante
Asume igual matriz de covarianza de las variables independientes entre cada grupo
Requiere de un gran tamaño muestral para obtener coeficientes insesgados (no ocurre así con la regresión logística)
Validación del LDA
Contraste de normalidad Shapiro-Wilk o kolmogorov-smirnov para cada variable
Contraste de normalidad multivariante de Henze-Zirkler, Doornik-Hansen o Royston con el paquete El paquete MVN
La homogeneidad de varianzas se puede realizar con el El test Box M (1949) del paquete biotools. El test Box M es muy sensible a violaciones de la normalidad multivariante, por lo que se recomienda emplear un límite de significancia de 0.001
Matriz de Confusión para comparar la clasificación estimada con la clasificación observada. El porcentaje de individuos correctamente clasificados nos va a dar información sobre la confianza de la ecuación predictiva.
Generalmente en Antropología Forense se considerará una ecuación confiable aquella que tenga al menos un 85% de casos correctamente clasificados.
EJEMPLO
Usaremos el mismo fichero craneos.xlsx con los datos del Tibet
library(readxl)
craneos <- read_excel("data/craneos.xlsx") #Ahora le llamamos craneos al objeto
#Revisamos la data
str(craneos)
tibble [32 x 6] (S3: tbl_df/tbl/data.frame) $ Longitud : num [1:32] 19050 17250 16700 16950 17500 ... $ Anchura : num [1:32] 15250 13200 13000 15050 13850 ... $ Altura : num [1:32] 14500 12550 12550 13350 12600 ... $ Altura_Cara : num [1:32] 7350 6300 6950 6450 7750 7150 7050 7350 7000 6200 ... $ Anchura_Cara: num [1:32] 13650 12100 11950 12800 13550 ... $ Tipo_raza : num [1:32] 1 1 1 1 1 1 1 1 1 1 ...
summary(craneos[ , c(1:5)])
Longitud Anchura Altura Altura_Cara Anchura_Cara Min. :16250 Min. :12650 Min. :12150 Min. :6200 Min. :11850 1st Qu.:17325 1st Qu.:13488 1st Qu.:12962 1st Qu.:6988 1st Qu.:13138 Median :17950 Median :13900 Median :13300 Median :7350 Median :13450 Mean :17994 Mean :13906 Mean :13330 Mean :7294 Mean :13370 3rd Qu.:18463 3rd Qu.:14250 3rd Qu.:13662 3rd Qu.:7650 3rd Qu.:13662 Max. :20000 Max. :15300 Max. :14500 Max. :8250 Max. :14650
Convertimos factor la variable Tipo_raza para asegurarnos que sea leida como categórica
craneos$Tipo_raza <- factor(craneos$Tipo_raza)
Realizamos los graficos QQ-plot para probar normalidad
par(mfrow=c(2,2))
for(i in 1:4){
qqnorm(unlist(craneos[,i]), pch = 1)
qqline(unlist(craneos[,i]), col = "steelblue", lwd = 2)}
par(mfrow=c(1,1))
Visualizamos el último gráfico
qqnorm(unlist(craneos[,5]), pch = 1) qqline(unlist(craneos[,5]), col = "steelblue", lwd = 2)
Ahora realizamos los Test de Normalidad para cada variable
by(craneos$Longitud, craneos$Tipo_raza, shapiro.test)
craneos$Tipo_raza: 1
Shapiro-Wilk normality test
data: dd[x, ]
W = 0.97442, p-value = 0.8901
------------------------------------------------------------
craneos$Tipo_raza: 2
Shapiro-Wilk normality test
data: dd[x, ]
W = 0.94141, p-value = 0.4004
by(craneos$Anchura, craneos$Tipo_raza, shapiro.test)
craneos$Tipo_raza: 1
Shapiro-Wilk normality test
data: dd[x, ]
W = 0.96081, p-value = 0.6468
------------------------------------------------------------
craneos$Tipo_raza: 2
Shapiro-Wilk normality test
data: dd[x, ]
W = 0.94619, p-value = 0.4666by(craneos$Altura, craneos$Tipo_raza, shapiro.test)
craneos$Tipo_raza: 1
Shapiro-Wilk normality test
data: dd[x, ]
W = 0.9586, p-value = 0.6053
------------------------------------------------------------
craneos$Tipo_raza: 2
Shapiro-Wilk normality test
data: dd[x, ]
W = 0.9654, p-value = 0.785
by(craneos$Altura_Cara, craneos$Tipo_raza, shapiro.test)
craneos$Tipo_raza: 1
Shapiro-Wilk normality test
data: dd[x, ]
W = 0.96785, p-value = 0.7798
------------------------------------------------------------
craneos$Tipo_raza: 2
Shapiro-Wilk normality test
data: dd[x, ]
W = 0.97646, p-value = 0.9396by(craneos$Anchura_Cara, craneos$Tipo_raza, shapiro.test)
craneos$Tipo_raza: 1
Shapiro-Wilk normality test
data: dd[x, ]
W = 0.94674, p-value = 0.4071
------------------------------------------------------------
craneos$Tipo_raza: 2
Shapiro-Wilk normality test
data: dd[x, ]
W = 0.93485, p-value = 0.322
Ahora realizamos el Test de Normalidad Multiple con el paquete MVN
library(MVN) mvn(craneos[ , c(1:5)],mvnTest = "hz") #Test de Henze-Zirkler
$multivariateNormality
Test HZ p value MVN
1 Henze-Zirkler 0.8261507 0.3755845 YES
$univariateNormality
Test Variable Statistic p value Normality
1 Anderson-Darling Longitud 0.4211 0.3048 YES
2 Anderson-Darling Anchura 0.3057 0.5470 YES
3 Anderson-Darling Altura 0.2250 0.8049 YES
4 Anderson-Darling Altura_Cara 0.1872 0.8962 YES
5 Anderson-Darling Anchura_Cara 0.7763 0.0392 NO
$Descriptives
n Mean Std.Dev Median Min Max 25th 75th Skew
Longitud 32 17993.75 936.5129 17950 16250 20000 17325.0 18462.5 0.4267349
Anchura 32 13906.25 684.1229 13900 12650 15300 13487.5 14250.0 0.2908175
Altura 32 13329.69 608.2576 13300 12150 14500 12962.5 13662.5 0.0930111
Altura_Cara 32 7293.75 539.0778 7350 6200 8250 6987.5 7650.0 -0.2261243
Anchura_Cara 32 13370.31 744.4265 13450 11850 14650 13137.5 13662.5 -0.4722848
Kurtosis
Longitud -0.6532942
Anchura -0.6738678
Altura -0.8202897
Altura_Cara -0.8175725
Anchura_Cara -0.3977022Aplicamos la función boxM() del paquete biotools que realiza el contraste de homogeneidad de matrices de covarianza.
#devtools::install_github("arsilva87/biotools") #Para instalar este paquete usamos este código
library(biotools)
boxM(data = craneos[, 1:5], grouping = craneos$Tipo_raza)
Box's M-test for Homogeneity of Covariance Matrices
data: craneos[, 1:5]
Chi-Sq (approx.) = 18.371, df = 15, p-value = 0.2437
Observamos que no tenemos elementos para rechazar la homogeneidad de covarianzas
Ahora que todos los test nos dan bien, procedemos a realizar el Analisis Discriminante Lineal con el paquete MASS
#install.packages("MASS") #Recuerden instalarlo primero
library(MASS)
modelo_lda <- lda(Tipo_raza ~ Longitud+Anchura+Altura+Altura_Cara+Anchura_Cara, data = craneos)
Podemos graficar el modelo
plot(modelo_lda)
Con los coefiencientes de la función discriminante, construimos la ecuación predictiva:
\(\hat Y=(0.0004772 * Longitud)+(-0.0008324*Anchura)+(-0.0000279*Altura)+\) \((0.0009469*\verb|Altura_Cara|)+(0.000948*\verb|Anchura_Cara|)\)
Solo nos falta la constante, que el paquete MASS no lo muestra automáticamente, aunque internamente si lo usa para los calculos. Nosotros podemos obtenerlo con el siguiente código:
groupmean<-(modelo_lda$prior %*% modelo_lda$means) #operador % * % es para realizar producto de dos matrices. constant<--(groupmean %*% modelo_lda$scaling) constant
LD1 [1,] -16.22159
Para ayudarnos a clasificar los valores o scores de la ecuación predictiva, usamos como referencia los centroides de las funciones. Para eso, sacamos los promedios de los scores predictivos
lda_prediction <- predict(modelo_lda) #Aqui obtenemos los scores y clases pronosticadas tapply(lda_prediction$x, craneos$Tipo_raza, mean) #Así obtenemos las funciones de los centroides de los grupos
1 2 -0.8771315 0.9940824
Los cálculos de los scores y clase pronosticada se realizan de la siguiente forma:
Todo esto se guarda en el objeto “lda_prediction”, y tiene 3 elementos importantes:
Para la exactitud del modelo, usaremos una tabla de contingencia o matriz de confusión para tabular la clasificación estimada con la clasificación observada y ver las coincidencias
table(predict=lda_prediction$class, Actual=craneos$Tipo_raza)
Actual
predict 1 2
1 14 3
2 3 12
Luego calculamos la Exactitud (Accuracy) como: (14+12/32)*100 = 81.25%.
mean(lda_prediction$class==craneos$Tipo_raza)*100
[1] 81.25
En resumen, el modelo clasificó correctamente el 81.25% de las observaciones
Graficos de Partición por cada par de variables
#install.packages("klaR") #Recuerden instalarlo primero
library(klaR)
partimat(craneos[, 1:5],grouping=as.factor(craneos$Tipo_raza), method = "lda")
Finalmente, graficamos, como en este caso tenemos una sola función discriminante, procedemos así:
prediction <- data.frame(lda_prediction) #Guardamos nuestro objeto de predicción como un data frame
ggplot(prediction, aes(x=as.factor(class) , y=LD1, fill = as.factor(class)))+
geom_dotplot(binaxis = "y", stackdir = "center", binwidth = 0.25)+
xlab("Pronostico") + theme_light() + theme(legend.position = "none")
EJEMPLO CON MÁS DE DOS GRUPOS
Para este ejemplo, tomaremos los datos de las frecuencias sanguíneas de ABO, disponibles en la página de Wikipedia
Posteriormente, se limpió y se agregó los continentes de cada país.
Blood_type_new <- read_excel("data/Blood_type_final.xlsx")
head(Blood_type_new)
# A tibble: 6 x 11 continent country population `O+` `A+` `B+` `AB+` `O-` `A-` `B-` <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> 1 Europa Albania 3074579 0.341 0.312 0.145 0.052 0.06 0.055 0.026 2 Africa Algeria 43576691 0.4 0.3 0.15 0.0425 0.066 0.023 0.011 3 America Argentina 45479118 0.489 0.315 0.08 0.0245 0.049 0.0316 0.008 4 Asia Armenia 3021324 0.29 0.463 0.12 0.056 0.02 0.037 0.01 5 Oceania Australia 25466459 0.38 0.32 0.12 0.04 0.07 0.06 0.02 6 Europa Austria 8859449 0.3 0.37 0.12 0.05 0.06 0.07 0.02 # i 1 more variable: `AB-` <dbl>
Nos aseguramos que la variable continente sea leida como categórica
Blood_type_new$continent <- factor(Blood_type_new$continent)
Y eliminamos a Oceanía por tener poca frecuencia
Blood_type_new <- Blood_type_new[Blood_type_new$continent != "Oceania", ]
Revisamos las variables para asegurarnos que no existan valores perdidos o extraños
summary(Blood_type_new[ , c(4:11)])
O+ A+ B+ AB+
Min. :0.2680 Min. :0.1400 Min. :0.0472 Min. :0.00500
1st Qu.:0.3247 1st Qu.:0.2479 1st Qu.:0.1000 1st Qu.:0.02900
Median :0.3817 Median :0.2990 Median :0.1510 Median :0.04295
Mean :0.4044 Mean :0.2956 Mean :0.1637 Mean :0.04756
3rd Qu.:0.4684 3rd Qu.:0.3507 3rd Qu.:0.2114 3rd Qu.:0.06313
Max. :0.7500 Max. :0.4630 Max. :0.3680 Max. :0.14400
O- A- B- AB-
Min. :0.00060 Min. :0.00040 Min. :0.00010 Min. :0.000100
1st Qu.:0.01798 1st Qu.:0.01000 1st Qu.:0.00585 1st Qu.:0.001225
Median :0.04000 Median :0.02715 Median :0.01275 Median :0.004150
Mean :0.03907 Mean :0.03284 Mean :0.01340 Mean :0.004937
3rd Qu.:0.06000 3rd Qu.:0.06000 3rd Qu.:0.02000 3rd Qu.:0.009025
Max. :0.13000 Max. :0.08000 Max. :0.03130 Max. :0.012000
Para trabajar menos variables, podemos sumar los grupos O, los grupos A, los grupos B y finalmente los AB, es decir, no tomar en cuenta los Rh.
Blood_type_new$O <- Blood_type_new$`O+`+Blood_type_new$`O−` Blood_type_new$A <- Blood_type_new$`A+`+Blood_type_new$`A−` Blood_type_new$B <- Blood_type_new$`B+`+Blood_type_new$`B−` Blood_type_new$AB <- Blood_type_new$`AB+`+Blood_type_new$`AB−`
Para efectos practicos vamos obviar las pruebas de normalidad. Pero quedan de tarea.
Realizamos el análisis discriminante
modelo_lda <- lda(continent ~ O+A+B+AB, data = Blood_type_new)
Warning in lda.default(x, grouping, ...): group Oceania is empty
modelo_lda
Call:
lda(continent ~ O + A + B + AB, data = Blood_type_new)
Prior probabilities of groups:
Africa America Asia Europa
0.1750000 0.1583333 0.3500000 0.3166667
Group means:
O A B AB
Africa 0.5032857 0.2672524 0.1932762 0.04118571
America 0.5735368 0.2982684 0.1036421 0.02454737
Asia 0.3974405 0.2963714 0.2395762 0.06745714
Europa 0.3963421 0.4127105 0.1356842 0.05617105
Coefficients of linear discriminants:
LD1 LD2 LD3
O 16.99326 30.00720 65.17075
A 14.39362 12.43382 74.58440
B 33.72184 19.74911 78.34340
AB 22.18339 21.79194 24.73584
Proportion of trace:
LD1 LD2 LD3
0.5183 0.4659 0.0158 Como tenemos más de dos funciones discriminantes, podemos graficar las dos primeras:
plot(modelo_lda, dimen=2, col = as.integer(Blood_type_new$continent))
Verificamos la exactitud del modelo
lda_prediction <- predict(modelo_lda) table(predict=lda_prediction$class, Actual=Blood_type_new$continent)
Actual predict Africa America Asia Europa Oceania Africa 8 2 5 0 0 America 3 13 0 2 0 Asia 8 0 28 3 0 Europa 2 4 9 33 0 Oceania 0 0 0 0 0
Entonces la Exactitud (Accuracy) es: (8+13+28+33/120)*100 = 68.3%%.
mean(lda_prediction$class==Blood_type_new$continent)*100
[1] 68.33333
El modelo clasifica correctamente el 68.3% de las observaciones
